package com.jidolstar.astro
{
	/**
	 * Precession(세차)<br>
	 * 세차운동에 관련된 변환행렬을 만들어준다.
	 * 적도좌표계와 황도좌표계에 대응한다. 
	 * IAU에서  1976년에 채택한 방법을 사용한다. 
	 * 
	 * @author Yongho,Ji (jidolstar@gmail.com)
	 * @since 2009.04.11
	 * 
	 * @see http://en.wikipedia.org/wiki/Precession_(astronomy)
	 * @see http://www.mps.mpg.de/homes/fraenz/systems/systems2art/node3.html
	 * @see http://articles.adsabs.harvard.edu/full/1979A%26A....73..282L
	 * @see Astronomical personal on the computer 4th p20~p24
	 * @see Astronomical Algorithms 2nd p134~138
	 * 
	 * @example
	 * <listing version="3.0">
	 *	//-----------------------------------------------------
	 *	// Astronomical Algorithms 2nd p.137
	 *	//-----------------------------------------------------			
	 *	var toJDE:Number = ALDate.date2JD( -214, 6, 30, 0, 0, 0 );
	 *	var fromJDE:Number = ALConst.JD_J2000;
	 *	var precession:ALPrecession = new ALPrecession( fromJDE, toJDE );
	 *	var eclMat:ALMatrix3D = precession.getPrecessionOfEclMatrix();
	 *	var eclLong:Number = 149.48194;
	 *	var eclLat:Number = 1.76549;
	 *	var eclVec:ALVector3D = new ALVector3D( eclLong * ALConst.D2R, eclLat * ALConst.D2R, 1, ALVector3D.COORDINATE_TYPE_SPHERICAL );
	 *	var precEclVec:ALVector3D =  eclMat.multiply3x3ByVector3D( [eclVec] )[0];
	 * 
	 *  trace( "------------------------------" );
	 *	trace( "λ0 = " + eclVec.longitudeDegree ); //λ0 = 149.48194
	 *	trace( "β0 = " + eclVec.latitudeDegree ); //β0 = 1.76549
	 *	trace( "fromJDE = " + fromJDE ); //fromJDE = 2451545
	 *	trace( "toJDE = " + toJDE ); //toJDE = 1643074.5
	 *	trace( "λ = " + precEclVec.longitudeDegree ); //λ = 118.70416774860918
	 * 	trace( "β = " + precEclVec.latitudeDegree ); //β = 1.6153320060500687
	 *	
	 *	//-----------------------------------------------------
	 *	// Astronomical Algorithms 2nd p.135
	 *	//-----------------------------------------------------
	 *	var RA:Number = ALHMS.HMS2FractionalHour( 2, 44, 11.986 ); //적경 2h 44m 11.986s
	 *	var DE:Number = ALDMS.DMS2FractionalDegree( 49,13,42.48 ); //적위 49°13'42".48
	 *	var pRA:Number = 0.03425;	//annual proper motion in right ascension [second]
	 *	var pDE:Number = -0.0895;	//annual proper motion in declination [arcsecond]
	 *	var fromDate:ALDate = new ALDate( ALConst.JD_J2000 );	//2000년 1월 1일 
	 *	var toDate:ALDate = new ALDate( 2028, 11, 13.19 );	//2028년 11월 13.19일
	 *	fromJDE = fromDate.JD;
	 *	toJDE = toDate.JD;
	 *	var julianYear:Number = ALTimeConverter.getJulianYear( fromJDE, toJDE ); //Julian Year
	 *	var deltaRA:Number = pRA * julianYear;
	 *	var deltaDE:Number = pDE * julianYear;
	 *	var RA1:Number = RA + deltaRA / 3600;
	 *	var DE1:Number = DE + deltaDE / 3600;
	 *	precession = new ALPrecession( fromJDE, toJDE );
	 *	var equVec:ALVector3D = new ALVector3D( RA1 * ALConst.H2R, DE1 * ALConst.D2R, 1, ALVector3D.COORDINATE_TYPE_SPHERICAL );
	 *	var equMat:ALMatrix3D = precession.getPrecessionOfEquMatrix();
	 *	var precEquVec:ALVector3D = equMat.multiply3x3ByVector3D( [equVec] )[0];
	 * 
	 *	trace( "------------------------------" );
	 *	trace( "α0 = " + RA * ALConst.H2D ); //α0 = 41.04994166666666
	 *	trace( "δ0 = " + DE); //δ0 = 49.22846666666667
	 *	trace( "fromJDE = " + fromJDE ); //fromJDE = 2451545
	 *	trace( "toJDE = " + toJDE ); //toJDE = 2462088.69
	 * 	trace( "Julian Year = " +julianYear); //28.86704996577671
	 *	trace( "Δα = " + deltaRA ); //Δα = 0.9886964613278525
	 *	trace( "Δδ = " + deltaDE ); //Δδ = -2.5836009719370154
	 *	trace( "α0 + Δα = " + equVec.longitudeDegree ); //α0 + Δα = 41.05406123525553
	 *	trace( "δ0 + Δδ = " + equVec.latitudeDegree ); //δ0 + Δδ = 49.22774899973002
	 *	trace( "α = " + precEquVec.longitudeDegree + " = " + ALHMS.fractionalHour2HMS( precEquVec.longitudeHour ).toString("h","m","s",3)); //α = 41.5472125975129 = 02h46m11.331s
	 *	trace( "δ = " + precEquVec.latitudeDegree + " = " + ALDMS.fractionalDegree2DMS( precEquVec.latitudeDegree ).toString());//δ = 49.34848211286583 = 49° 20’ 54”.54	  
	 * </listing>
  	 */ 
	public class ALPrecession
	{
		static public const BASED_SYSTEM_IAU1976:String = "IAU1976";
		static public const BASED_SYSTEM_IAU2000A:String = "IAU2000A";
		static public const BASED_SYSTEM_IAU2006:String ="IAU2006";
		private var t:Number;
		private var T:Number;
		private var _fromJDE:Number;
		private var _toJDE:Number;
		private var _basedSystem:String;
		
		/**
		 * 생성자
		 * @param fromJDE Epoch given 
		 * @param toJDE Epoch to precess to
		 */ 
		public function ALPrecession( fromJDE:Number, toJDE:Number, basedSystem:String = "IAU1976" )
		{
			setJDE( fromJDE, toJDE );
		}
		
		/**
		 * Julian Ephemeris Date(JDE) 설정 
		 * @param fromJDE Epoch given 
		 * @param toJDE Epoch to precess to
		 */ 
		public function setJDE( fromJDE:Number, toJDE:Number, basedSystem:String = "IAU1976" ):void
		{
			switch( basedSystem )
			{
				case BASED_SYSTEM_IAU1976:
				case BASED_SYSTEM_IAU2000A:
					break;
				default:
					basedSystem = BASED_SYSTEM_IAU1976; 
					break;
			}
			_fromJDE = fromJDE;
			_toJDE = toJDE;
			_basedSystem = basedSystem;
			
			T = ALTimeConverter.getJulianCenturyJ2000( fromJDE );
			t = ALTimeConverter.getJulianCentury( fromJDE, toJDE );
		}
		
		/**
		 * 세차운동시 사용하는 기준 시스템. 기본값은 IAU1976 이다.
		 */ 
		public function get basedSystem():String
		{
			return _basedSystem;
		}
		
		
		/**
		 * Epoch 
		 */ 
		public function get fromJDE():Number
		{
			return _fromJDE;
		}
	
		/**
		 * Epoch to precess to 
		 */
		public function get toJDE():Number
		{
			return _toJDE;
		}
		
		/**
		 * Precession of Equatorial Coordinate Transformation Matrix
		 * @return 세차운동이 적용된 적도좌표변환 행렬  
		 */ 
		public function getPrecessionOfEquMatrix():ALMatrix3D
		{
			var zeta:Number; //ζ
			var z:Number; //z 
			var theta:Number; ////ϑ
			
			switch( basedSystem )
			{
				case BASED_SYSTEM_IAU1976:
					zeta = ( (2306.2181+(1.39656-0.000139*T)*T)+
			               ((0.30188-0.000344*T)+0.017998*t)*t )*t* ALConst.S2R;
			        z = zeta + ( (0.79280+0.000411*T)+0.000205*t)*t*t* ALConst.S2R;
			        //z = ( (2306.2181+(1.39656-0.000139*T)*T)+((1.09468+0.000066*T)+0.018203*t)*t )*t* ALConst.S2R;
			        theta = ( (2004.3109-(0.85330+0.000217*T)*T)-
			                ((0.42665+0.000217*T)+0.041833*t)*t )*t* ALConst.S2R;
					break;
				case BASED_SYSTEM_IAU2000A:
					break;
				case BASED_SYSTEM_IAU2006:
					break;
			}

			
			//trace( " ζ = " + ALMath.mapTo0To2PIRange( zeta ) * ALConst.R2D); 
			//trace( " z = " + ALMath.mapTo0To2PIRange( z ) * ALConst.R2D );
			//trace( " ϑ = " + ALMath.mapTo0To2PIRange( theta ) * ALConst.R2D);

			/*
			이 식은 아래 식과 동일하다.
			var A:ALMatrix3D = new ALMatrix3D();
			var B:ALMatrix3D = new ALMatrix3D();
			var C:ALMatrix3D = new ALMatrix3D();
			A.setZAxisRotation( -z );
			B.setYAxisRotation( theta );
			C.setZAxisRotation( -zeta );
			return ALMatrix3D.multiplyMatrix( ALMatrix3D.multiplyMatrix( A, B ), C );  
			*/	
			
			var cos_zeta:Number = Math.cos( zeta );
			var sin_zeta:Number = Math.sin( zeta );
			var cos_z:Number = Math.cos( z );
			var sin_z:Number = Math.sin( z );
			var cos_t:Number = Math.cos( theta );
			var sin_t:Number = Math.sin( theta );
			return new ALMatrix3D( 	cos_z * cos_t * cos_zeta - sin_z * sin_zeta,
									-cos_z * cos_t * sin_zeta - cos_zeta * sin_z,
									-cos_z * sin_t,
									sin_z * cos_t * cos_zeta + cos_z * sin_zeta,
									-sin_z * cos_t * sin_zeta + cos_z * cos_zeta,
									-sin_t * sin_z,
									sin_t * cos_zeta,
									-sin_t * sin_zeta,
									cos_t ); 				  			
		}

		/**
		 * Precession of Ecliptic Coordinate Transformation Matrix
		 * @return 세차운동이 적용된 황도좌표변환 행렬  
		 */ 
		public function getPrecessionOfEclMatrix():ALMatrix3D
		{
			var P:Number; //Π : ascending node longitude of the plane of date  relative to the plane of date 
			var p:Number; //π : The orientation of the ecliptic plane of date with respect to the the ecliptic plane at another date is defined by the inclination 
			var eta:Number; //η : the angle between the ecliptic at the starting epoch and the ecliptic at the final epoch
			switch( basedSystem )
			{
				case BASED_SYSTEM_IAU1976:
					P  = 174.876383889* ALConst.D2R + 
		        		( ((3289.4789+0.60622*T)*T) +
		              	((-869.8089-0.50491*T) + 0.03536*t)*t ) * ALConst.S2R;
					p  = ( (5029.0966+(2.22226-0.000042*T)*T)+
					   	((1.11113-0.000042*T)-0.000006*t)*t )*t * ALConst.S2R;
					eta = ( (47.0029-(0.06603-0.000598*T)*T)+
					   	( (-0.03302+0.000598*T)+0.000060*t)*t )*t * ALConst.S2R;
					break;
				case BASED_SYSTEM_IAU2000A:
					break;
				case BASED_SYSTEM_IAU2006:
					break;
			}			
			
			//trace( " η = " + (ALMath.mapTo0To2PIRange( eta ) * ALConst.R2D - 360)); 
			//trace( " Π = " + ALMath.mapTo0To2PIRange( P ) * ALConst.R2D );
			//trace( " π = " + (ALMath.mapTo0To2PIRange( p ) * ALConst.R2D - 360));
			
			
			/*
			이 식은 아래 식과 동일하다.
			var A:ALMatrix3D = new ALMatrix3D();
			var B:ALMatrix3D = new ALMatrix3D();
			var C:ALMatrix3D = new ALMatrix3D();
			A.setZAxisRotation( -(P+p) );
			B.setXAxisRotation( eta );
			C.setZAxisRotation( P );
			return ALMatrix3D.multiplyMatrix( ALMatrix3D.multiplyMatrix( A, B ), C );  
			*/
			
			var ca:Number = Math.cos( -P-p );
			var sa:Number = Math.sin( -P-p );	
			var cb:Number = Math.cos( eta );
			var sb:Number = Math.sin( eta );	
			var cc:Number = Math.cos( P );
			var sc:Number = Math.sin( P );
			
			return new ALMatrix3D(  ca * cc - sa * cb * sc,
									ca * sc + sa * cb * cc,
									sa * sb,
									-sa * cc - ca * cb * sc,
									-sa * sc + ca * cb * cc,
									ca * sb,
									sb * sc,
									-sb * cc,
									cb );				
		}

	}
}